2019-02-11

Temario

  1. Introducción
  2. Regresión Lineal
    • Simple
    • Múltiple
  3. Estudios observacionales
    • Transversales
    • Caso-Control
  4. Modelos Lineales Generalizados
    • Odds Ratio
    • Prevalence Ratio
  5. Estudios adicionales
    • Cohortes - Risk Ratio
    • Tiempo a evento - Hazard Ratio
    • Desenlace conteo o tasa - Incidence Risk Ratio

Introducción

¿qué es R?

  • Software Libre

¿por qué usar R?

  • Ciencia Reproducible:
    • limpiar bases de datos
    • visualizar
    • ejecutar modelos
    • comunicar resultados

¿cómo usar R?

  • IDE: RStudio
  • Paquetes
  • Funciones

Introducción

¿por qué usar una regresión?

  • permite fijar al menos una variable independiente o exposición y observar una respuesta en la variable dependiente o desenlace.

  • permite explicar el cambio promedio de un evento \(Y\) en base a cambios en \(X\), usando coeficientes o medidas de asociación.

  • permite predecir la probabilidad asociada a un evento.

  • .

Regresión Lineal Simple

características

  • una variable independiente (simple)
  • una variable dependiente (univariada)
  • ambas variables deben ser numéricas

objetivos

  • ajustar datos a una recta
  • interpretar medida de bondad de ajuste ( R2 ) y coeficientes
  • evaluar supuestos
  • visualizar el modelo

\[ Y = \beta_0 + \beta_1 X_1 + \epsilon \]

RLS: datos

espir <- read_dta("data-raw/espirometria.dta") %>% as_factor()
espir
## # A tibble: 654 x 7
##     caso codigo  edad   vef talla sexo   fumar
##    <dbl>  <dbl> <dbl> <dbl> <dbl> <fct>  <fct>
##  1     1    301     9  1.71  57   female No   
##  2     2    451     8  1.72  67.5 female No   
##  3     3    501     7  1.72  54.5 female No   
##  4     4    642     9  1.56  53   male   No   
##  5     5    901     9  1.89  57   male   No   
##  6     6   1701     8  2.34  61   female No   
##  7     7   1752     6  1.92  58   female No   
##  8     8   1753     6  1.41  56   female No   
##  9     9   1901     8  1.99  58.5 female No   
## 10    10   1951     9  1.94  60   female No   
## # ... with 644 more rows

RLS: distribución

  • supuestos:
    • #1 linealidad
    • #2 independencia de observaciones
espir %>% 
  ggplot(aes(edad,vef)) +
  geom_point()

  • error = residual

RLS: suma de mínimos cuadrados

Cálculo de la sumatoria del cuadrado de los residuales hacia la media y la recta:

\[ SSE(mean) = \sum (data - mean)^2 \] \[ SSE(fit) = \sum (data - fit)^2 \]

\[ Var(x) = \frac{SSE(x)^2}{n} \]

Medida de bondad de ajuste:

\[ R^2 = \frac{Var(mean) - Var(fit)}{Var(mean)} \]

RLS: R2

wm1 <- lm(vef ~ edad, data = espir)
wm1 %>% glance() %>% select(r.squared:df)
## # A tibble: 1 x 6
##   r.squared adj.r.squared sigma statistic   p.value    df
## *     <dbl>         <dbl> <dbl>     <dbl>     <dbl> <int>
## 1     0.572         0.572 0.568      872. 2.45e-122     2

INTERPRETACIÓN

  • edad explica el 57% de la variabilidad de VEF
  • existe un 57% de reducción en la variabilidad de VEF al tomar en cuenta la edad

RLS: coeficientes

wm1 %>% tidy()
## # A tibble: 2 x 5
##   term        estimate std.error statistic   p.value
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>
## 1 (Intercept)    0.432   0.0779       5.54 4.36e-  8
## 2 edad           0.222   0.00752     29.5  2.45e-122
wm1 %>% confint_tidy()
## # A tibble: 2 x 2
##   conf.low conf.high
##      <dbl>     <dbl>
## 1    0.279     0.585
## 2    0.207     0.237
  • \(\beta_{edad}\): En la población,
  • por cada incremento de edad en una unidad, el VEF en promedio incrementa en 0.22 mL/s,
  • con un intervalo de confianza AL 95% de 0.21 a 0.24 mL/s.
  • Este resultado es estadísticamente significativo con un valor p < 0.001

RLS: supuesto #3 normalidad

RLS: supuesto #3 normalidad

wm1 %>% 
  augment() %>% 
  ggplot() +
  geom_qq(aes(sample=.std.resid)) + 
  geom_qq_line(aes(sample=.std.resid))

  • META: todos los puntos sobre la línea
wm1 %>% augment() %>% 
  ggplot(aes(.std.resid)) +
  geom_histogram()

RLS: supuesto #4 homoscedasticidad

RLS: supuesto #4 homoscedasticidad

  • META: distribución idéntica a ambos lados de la línea
wm1 %>% 
  augment() %>% 
  ggplot(aes(.fitted,.std.resid)) +
  geom_point() +
  geom_smooth() +
  geom_hline(yintercept = c(0))

RLS: ¿cómo se ve el modelo?

\[ Y = \beta_0 + \beta_1 X_1 + \epsilon \] \[ VEF = 0.60 + 0.22 (edad) + \epsilon \]

espir %>% 
  ggplot(aes(edad,vef)) +
  geom_point() + 
  geom_smooth(method = "lm")

RLS: retroalimentación

  • R2 indica el porcentaje de variabilidad del desenlace (var. dependiente) explicada por el predictor (var. independiente).

  • los coeficientes permiten predecir el posible valor del desenlace en base a un modelo estadístico.

  • los supuestos permiten evaluar qué tan adecuado es el ajuste de los datos al modelo.

RLS: ejercicio

  • ajustar una recta
espir %>% 
  ggplot(aes(edad,talla)) +
  geom____()
  • identificar coeficientes y R2
# recordar: y ~ x
wm1 <- lm(_____ ~ _____, data = _____)
wm1 %>% g_____
wm1 %>% t_____
wm1 %>% c_____
  • evaluar supuestos: normalidad y homoscedasticidad
wm1 %>% augment() %>% 
  ggplot() + _____

wm1 %>% augment() %>% 
  ggplot(aes(.fitted,.std.resid)) + __________

Intermedio #1: Tabla 1 y 2

“… ¡en un solo comando!”

compareGroups(fumar ~ edad + vef + talla + sexo,
              data = espir, byrow=T#,method=c(vef=2)
              ) %>% 
  createTable(show.all = T) #%>% export2xls("table/tab1.xls")
## 
## --------Summary descriptives table by 'Habitos de fumar'---------
## 
## ________________________________________________________ 
##               [ALL]        No          Si      p.overall 
##               N=654       N=589       N=65               
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯ 
## edad       9.93 (2.95) 9.53 (2.74) 13.5 (2.34)  <0.001   
## vef        2.64 (0.87) 2.57 (0.85) 3.28 (0.75)  <0.001   
## talla      61.1 (5.70) 60.6 (5.67) 66.0 (3.19)  <0.001   
## sexo:                                            0.071   
##     male   336 (51.4%) 310 (92.3%) 26 (7.74%)            
##     female 318 (48.6%) 279 (87.7%) 39 (12.3%)            
## ¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯

Intermedio: distribuciones

espir %>% 
  select(edad,vef,talla) %>% 
  gather(key,value) %>% 
  ggplot(aes(value)) +
  geom_histogram() +
  facet_wrap(~key,scales = "free")

## # A tibble: 1 x 10
##   n_obs   min   max  mean    sd   q25   q50   q75 skewness kurtosis
##   <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
## 1   654 0.791  5.79  2.64 0.867  1.98  2.55  3.12    0.661     3.30

Regresión Lineal Múltiple

características

  • dos o más variables independientes (múltiple)
  • una variable dependiente numérica

objetivo

  • epidemiológico: controlar confusión en el análisis.
    • estandarización (directa/indirecta),
    • estratificación,
    • puntajes de propensión,
    • regresión multivariable
  • estadístico: emplear múltiples predictores
    • lineales, no-lineales, segmentadas

RLM: distribución

espir %>% 
  ggplot(aes(edad,vef,colour=sexo)) +
  geom_point()

RLM: R2 y coeficientes

## # A tibble: 1 x 6
##   r.squared adj.r.squared sigma statistic   p.value    df
## *     <dbl>         <dbl> <dbl>     <dbl>     <dbl> <int>
## 1     0.607         0.606 0.544      503. 9.49e-133     3
  • edad y sexo explica el 60% de la variabilidad de VEF.
## # A tibble: 3 x 7
##   term        estimate std.error statistic   p.value conf.low conf.high
##   <chr>          <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
## 1 (Intercept)    0.605   0.0781       7.74 3.79e- 14    0.451     0.758
## 2 edad           0.220   0.00722     30.6  7.32e-128    0.206     0.235
## 3 sexofemale    -0.323   0.0426      -7.59 1.13e- 13   -0.407    -0.240
  • \(\beta_{sex:female}\)
  • En la población, el VEF en promedio en mujeres es 0.32 mL/s menor que en hombres,
  • con un intervalo de confianza AL 95% de 0.24 a 0.41 mL/s, ajustando por edad.
  • Este resultado es estadísticamente significativo con un valor p < 0.001

RLM: supuesto normalidad y homoscedasticidad

wm1 %>% 
  augment() %>% 
  ggplot() +
  geom_qq(aes(sample=.std.resid)) + 
  geom_qq_line(aes(sample=.std.resid))

wm1 %>% 
  augment() %>% 
  ggplot(aes(.fitted,.std.resid)) +
  geom_point() +
  geom_smooth() +
  geom_hline(yintercept = c(0))

RLM: post-estimación

  • multicolinealidad: VIF

  • outliers: puntos influenciales o ruido

RLM: ¿cómo se ve el modelo?

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \epsilon \] \[ VEF = 0.60 + 0.22 (edad) - 0.32 (sexo) + \epsilon \]

wm1 %>% augment() %>% 
  ggplot(aes(colour= sexo)) +
  geom_point(aes(edad,vef),alpha=0.05) +
  geom_line(aes(edad,.fitted), lwd=1.5)

RLM: ejercicio

  • 3 predictores
wm1 <- lm(vef ~ edad + sexo + fumar, data = espir)
  • R2 y coeficientes
wm1 %>% g___()
wm1 %>% t___() 
wm1 %>% c___()
  • supuesto normalidad y homoscedasticidad
wm1 %>% augment() %>% 
  ggplot() +
  geom_qq________ + 
  geom________

wm1 %>% augment() %>% 
  ggplot(aes(_____,_____)) +
  geom_p____()

Intermedio #2: Modificación de efecto

¿cómo se ve el modelo?

\[ Y = \beta_0 + \beta_1 X_1 + \beta_2 X_2 + \beta_3 (X_1 * X_2) + \epsilon \]

## # A tibble: 4 x 4
##   term         estimate conf.low conf.high
##   <chr>           <dbl>    <dbl>     <dbl>
## 1 (Intercept)     0.253   0.0911     0.416
## 2 edad            0.243   0.226      0.259
## 3 fumarSi         1.94    1.13       2.76 
## 4 edad:fumarSi   -0.163  -0.223     -0.102
  • \(\beta_{interacción}\)
  • La diferencia en la media de VEF asociado al incremento en 1 año de edad en fumadores
  • menos
  • la diferencia en la media de VEF asociado al incremento en 1 año de edad en no fumadores
  • es -0.16 mL/s, con un IC 95% de -0.22 a -0.1.

Estudios Observacionales

  • Asociar:
    • condición (expuesto o no-expuesto) y
    • resultado (caso o control)
  • El tipo de análisis depende del tipo de estudio:
    • “Reclutamos \(x_1\) sujetos expuestos y \(x_2\) sujetos no expuestos
    • Tipo Cohorte

    • “Reclutamos \(x_1\) sujetos caso y \(x_2\) sujetos control
    • Tipo Caso-Control

tipo de estudio medida de asociación
caso - control OR
cohorte RR
transversal PR
tiempo a evento HR

Bases de datos

  • Cohorte
riskdb <- read_dta("data-raw/epidb.dta") %>% as_factor() %>% 
  mutate(aids=if_else(aids==0,"No","Yes") %>% as.factor(),
         ccr2=if_else(ccr2==1,"No","Yes") %>% as.factor())
  • Encuesta a nivel nacional
logdb <- read_dta("data-raw/brfss2000.dta") %>% as_factor() %>% 
  filter(!is.na(depress))
  #mutate(dn_cro=as.factor(dn_cro))

Cohorte

  • Elegimos los totales de expuestos y no expuestos
riskdb %>% 
  count(ccr2,aids) %>% 
  spread(aids,n) %>% 
  mutate(risk=Yes/(No+Yes)) %>% 
  mutate(rr=.$risk[2]/.$risk[1])
## # A tibble: 2 x 5
##   ccr2     No   Yes  risk    rr
##   <fct> <int> <int> <dbl> <dbl>
## 1 No      222   149 0.402 0.821
## 2 Yes      63    31 0.330 0.821
  • Si cambiamos el tamaño del grupo (no-expuestos x 100), no cambia el Riesgo Relativo (RR)
## # A tibble: 2 x 5
##   ccr2     No   Yes  risk    rr
##   <fct> <dbl> <dbl> <dbl> <dbl>
## 1 No    22200 14900 0.402 0.821
## 2 Yes      63    31 0.330 0.821

Caso Control

  • Elegimos los totales de casos y controles

  • Si cambiamos el tamaño del grupo (casos x 100), sí cambia el Riesgo Relativo (RR)

riskdb %>% 
  count(ccr2,aids) %>% 
  spread(aids,n) %>% 
  mutate(Yes=Yes*100,
         risk=Yes/(No+Yes)) %>% 
  mutate(rr=.$risk[2]/.$risk[1])
  • En estudios caso-control usamos el Odds Ratio (OR)
## # A tibble: 2 x 7
##   ccr2     No   Yes  risk  odds    rr    or
##   <fct> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 No      222 14900 0.985  67.1 0.995 0.733
## 2 Yes      63  3100 0.980  49.2 0.995 0.733
  • sobreestima el Riesgo Relativo (RR)

Modelos Lineales Generalizados

careacterística

objetivo

GLM: familias distribución y funciones de enlace

  • tabla

GLM: verosimilitud

GLM: Regresión Logística

  • Se aplica cuando el desenlace es categórico dicotómico

  • El término logístico refiere a la función de enlace logit, que convierte el resultado binario en una variable contínua de log(odds)

\[ log\left(\frac{\pi}{1-\pi}\right) = \beta_0 + \beta_1 X_1 + . . . + \beta_n X_n + \epsilon \]

  • El valor exponenciado de los coeficientes se pueden interpretar como Odds Ratio (OR)
data_frame(p=seq(0,1,0.01)) %>% 
  mutate(odd=p/(1-p)) %>% 
  ggplot(aes(p,log(odd))) +
  geom_point() +
  coord_flip()

GLM: ¿bondad de ajuste?

  • LogLikelihood
  • Information Criterias
    • AIC
    • BIC

GLM: ¿Supuestos?

OR: Caso-Control

  • Simularemos un escenario Caso-Control
  • Seleccionaremos al azar 500 casos y 500 controles
riskdb_or <- riskdb %>% 
  group_by(aids) %>% 
  sample_n(50) %>% 
  ungroup()
riskdb_or %>% ggplot(aes(aids,fill=ccr2)) +
  geom_bar(position = "fill")

riskdb_or %>% 
  count(aids) %>% 
  mutate(prc=n/sum(n))
## # A tibble: 2 x 3
##   aids      n   prc
##   <fct> <int> <dbl>
## 1 No       50   0.5
## 2 Yes      50   0.5
## # A tibble: 2 x 5
##   ccr2     No   Yes  odds    or
##   <fct> <int> <int> <dbl> <dbl>
## 1 No       41    45 1.10  0.506
## 2 Yes       9     5 0.556 0.506

OR: coeficientes de variables categóricas

wm1 <- glm(aids ~ ccr2, data = riskdb_or, family = binomial(link = "logit"))
## # A tibble: 2 x 6
##   term        log.odds  odds conf.low conf.high p.value
##   <chr>          <dbl> <dbl>    <dbl>     <dbl>   <dbl>
## 1 (Intercept)   0.0931 1.10     0.719      1.68   0.666
## 2 ccr2Yes      -0.681  0.506    0.145      1.59   0.255
  • \(\beta_0\)
  • En la población, el odds de tener sida dado que no poseen CCR2 es 1.1,
  • con un intervalo de confianza al 95% de 0.72 a 1.68

  • \(\beta_{CCR2:Yes}\)
  • En la población, el odds de tener sida dado que sí poseen CCR2 es
  • 0.51 veces,
  • el odds de tener depresión dado que no poseen CCR2,
  • con un intervalo de confianza al 95% de 0.15 a 1.59.
  • Este resultado no es estadísticamente significativo con un valor p = 0.255

OR: coeficientes de variables contínuas

wm1 <- glm(aids ~ age, data = riskdb_or, family = binomial(link = "logit"))
## # A tibble: 2 x 6
##   term        log.odds  odds conf.low conf.high p.value
##   <chr>          <dbl> <dbl>    <dbl>     <dbl>   <dbl>
## 1 (Intercept)   1.58   4.83     0.797     32.3   0.0929
## 2 age          -0.0451 0.956    0.906      1.01  0.0860
  • \(\beta_{edad}\)
  • En la población, por cada incremento en un año de edad,
  • el odds de sida cambia en 0.96, con un intervalo de confianza AL 95% de 0.91 a 1.01.
  • Este resultado no es estadísticamente significativo con un valor p = 0.086

RR: Cohorte

wm1 <- glm(aids ~ ccr2, data = riskdb, family = binomial(link = "log"))
## # A tibble: 2 x 6
##   term        log.odds  odds conf.low conf.high  p.value
##   <chr>          <dbl> <dbl>    <dbl>     <dbl>    <dbl>
## 1 (Intercept)   -0.912 0.402    0.353     0.452 5.48e-47
## 2 ccr2Yes       -0.197 0.821    0.586     1.10  2.18e- 1
  • \(\beta_0\)
  • En la población, el riesgo de sida dado que no poseen CCR2 es 0.4,
  • con un intervalo de confianza AL 95% de 0.35 a 0.45

  • \(\beta_{ccr2:yes}\)
  • En la población, el riesgo de tener sida dado que sí poseen CCR2 es
  • 0.82 veces,
  • el riesgo de tener sida dado que no poseen CCR2,
  • con un intervalo de confianza AL 95% de 0.59 a 1.1.
  • Este resultado no es estadísticamente significativo con un valor p = 0.218

PR: Transversal

wm1 <- glm(depress ~ short, data = logdb, family = binomial(link = "log"))
## # A tibble: 2 x 6
##   term        log.odds  odds conf.low conf.high p.value
##   <chr>          <dbl> <dbl>    <dbl>     <dbl>   <dbl>
## 1 (Intercept)   -2.23  0.107    0.103     0.112 0      
## 2 shortYes       0.121 1.13     1.04      1.23  0.00526
  • \(\beta_0\)
  • En la población, la prevalencia de depresión en sujetos que estatura alta es 0.11,
  • con un intervalo de confianza AL 95% de 0.1 a 0.11

  • \(\beta_{estatura:baja}\)
  • En la población, la prevalencia de tener depresión dado que se es de estatura baja es
  • 1.13 veces,
  • el prevalencia de tener depresión dado que se es de estatura alta,
  • con un intervalo de confianza AL 95% de 1.04 a 1.23.
  • Este resultado es estadísticamente significativo con un valor p = 0.005

HR: Tiempo a evento

Log-Rank test

HR: Tiempo a evento

Regresión de Cox

res.cox <- coxph(Surv(ttoaidsorexit, aids) ~ ccr2, data =  coxdb)
## # A tibble: 1 x 7
##   term    estimate std.error statistic p.value conf.low conf.high
##   <chr>      <dbl>     <dbl>     <dbl>   <dbl>    <dbl>     <dbl>
## 1 ccr2Yes   -0.327     0.198     -1.66  0.0978   -0.714    0.0602

IRR: Desenlaces conteo o tasas

glm(response~predictor1+predictor2+predictor3+ ... + offset(log(population)),
     family=poisson,data=...)
#https://stats.stackexchange.com/questions/66791/where-does-the-offset-go-in-poisson-negative-binomial-regression/66878#66878
#https://stackoverflow.com/questions/16046726/regression-for-a-rate-variable-in-r

Referencias

Presentaciones

  • Curso de Bioestadística. Maestría en Ciencias de la Investigación Epidemiológica. UPCH 2018.
  • José Alfredo Zavala. Regresión Lineal. 2018
  • Harold Akehurst. Bioestadística con R. Julio 2016

enlaces web

  • StatQuest: Linear Models Pt.1 - Linear Regression. Link: https://youtu.be/nk2CQITm_eo
  • Statistics for Biologist: Points of Significance. Link: https://www.nature.com/collections/qghhqm/pointsofsignificance